/*
 *  meshmanager.cpp
 *  gmit
 *
 *  Created by DONG Li on 4/11/10.
 *  Copyright 2010 LASG/IAP. All rights reserved.
 *
 */

#include <stdexcept>

#include "meshmanager.h"

template <typename VT, typename ET, typename CT>
MeshManager<VT, ET, CT>::MeshManager() {}

template <typename VT, typename ET, typename CT>
MeshManager<VT, ET, CT>::~MeshManager() {}

// **************************************************************** mesh file IO

template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::readMeshFile(const string& meshFileName) {}

template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::writeMeshFile() const
{
    cout << endl << "number of vertex grids: " << this->vertexGrid.size() << endl;
    for (int i = 0; i < this->vertexGrid.size(); ++i) {
        Grid<VT> vertex = this->vertexGrid[i];
        cout << "vertex " << i << ": (";
        const vector<double> x = vertex.getLocation();
        for (int j = 0; j < this->numDim-1; ++j) {
            cout << x[j] << ", ";
        }
        cout << x[this->numDim-1] << ")" << endl;
    }
    cout << endl << "number of vertex cells: " << this->vertexGridMesh.size() << endl;
    for (int i = 0; i < this->vertexGridMesh.size(); ++i) {
        Cell cell = this->vertexGridMesh[i];
        cout << "vertex cell " << i << ":" << endl;
        for (int j = 0; j < cell.getNumEdge(); ++j) {
            int edgeId = cell.getEdgeGlobalId(j);
            Edge edge = this->vertexGridMeshEdges[edgeId];
            cout << "  edge " << j << ": global id - " << edgeId << endl;
            for (int k = 0; k < edge.getNumVertex(); ++k) {
                int vertexId = edge.getVertexGlobalId(k);
                Grid<VT> vertex = this->vertexGrid[vertexId];
                const vector<double> x = vertex.getLocation();
                cout << "    vertex " << k << ": global id - " << vertexId;
                cout << " coordinate - (";
                for (int l = 0; l < this->numDim-1; ++l) {
                    cout << x[l] << ", ";
                }
                cout << x[this->numDim-1] << ")" << endl;
            }
        }
    }
}

// *************************************************** mesh construction methods
#ifdef MESHCREATOR
template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::setNumDim(const int nd)
{
    if (nd > 3 || nd < 1)
        throw invalid_argument("MeshManager::setNumDim: "
                               "Dimension must be 1, 2, or 3.");
    this->numDim = nd;
}

template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::setNumCell(const int dimId, const int n)
{
    if (this->meshStructure != Structured || this->meshForm != Uniform)
        throw runtime_error("MeshManager::setNumCell: "
                            "Invalid for nstructured mesh.");
    if (dimId < 0)
        throw invalid_argument("MeshManager::setNumCell: "
                               "Dimension id must not be less than 0.");
    if (dimId >= this->numDim)
        throw runtime_error("MeshManager::setNumCell: Wrong dimension id.");
    this->numCell[dimId] = n;
}

template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::setBndX(const int dimId,
                                      const double x1, const double x2)
{
    if (this->meshStructure != Structured || this->meshForm != Uniform)
        throw invalid_argument("MeshManager::setBndX: "
                               "Invalid for unstructured mesh.");
    if (dimId < 0)
        throw invalid_argument("MeshManager::setBndX: "
                               "Dimension id must not be less than 0.");
    if (dimId >= this->numDim)
        throw invalid_argument("MeshManager::setBndX: Wrong dimension id.");
    if (x1 >= x2 || x1 < 0)
        throw invalid_argument("MeshManager::setBndX: Wrong arguments.");
    this->bndX[dimId][0] = x1;
    this->bndX[dimId][1] = x2;
    this->span[dimId] = x2-x1;
    this->dx[dimId] = this->span[dimId]/this->numCell[dimId];
}

template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::setBndType(const int dimId, 
                                         const int sideId, const BoundaryType b)
{
    if (this->meshStructure != Structured || this->meshForm != Uniform)
        throw invalid_argument("MeshManager::setBndType: "
                               "Invalid for unstructured mesh.");
    if (dimId < 0)
        throw invalid_argument("MeshManager::setBndType: "
                               "Dimension id must not be less than 0.");
    if (dimId >= this->numDim)
        throw invalid_argument("MeshManager::setBndType: Wrong dimension id.");
    if (sideId < 0 || sideId > 1)
        throw invalid_argument("MeshManager::setBndType: Wrong side id.");
    this->bndType[dimId][sideId] = b;
}

template <typename VT, typename ET, typename CT>
double MeshManager<VT, ET, CT>::getDx(const int dimId) const
{
    if (this->meshStructure != Structured || this->meshForm != Uniform)
        throw runtime_error("getDx should not be called for unstructured mesh");
    if (dimId < 0)
        throw invalid_argument("Dimension id must not be less than 0.");
    if (dimId >= this->numDim)
        throw runtime_error("Wrong call of getDx, check dimension.");
    return this->dx[dimId];
}

template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::
createVertexGridCell(const CellType cellType,
                     const vector<Grid<VT> >& vertices,
                     int* vertexConn)
{
    Cell cell;
    // set vertices
    for (int i = 0; i < vertices.size(); ++i) {
        // check vertex's existence
        int vertexGlobalId;
        bool nonexist = true;
        for (int j = 0; j < this->vertexGrid.size(); ++j) // This should be optimized
            if (vertices[i].isDuplicate(this->vertexGrid[j])) {
                nonexist = false; // vertex exists, use it
                vertexGlobalId = j;
                break;
            }
        if (nonexist) { // vertex doesn't exsit, add it
            vertexGlobalId = this->vertexGrid.size();
            this->vertexGrid.push_back(vertices[i]);
        }
        cell.addVertexGlobalId(vertexGlobalId);
    }
    // set edges
    int numEdge;
    switch (cellType) {
        case Triangular:
            numEdge = 3;
            break;
        case Quadrilateral:
            numEdge = 4;
            break;
        default:
            break;
    }
    for (int i = 0; i < numEdge; ++i) {
        int localId1, localId2;
        int globalId1, globalId2;
        localId1 = *(vertexConn+i*2);
        localId2 = *(vertexConn+i*2+1);
        globalId1 = cell.getVertexGlobalId(localId1);
        globalId2 = cell.getVertexGlobalId(localId2);
        Edge edge;
        edge.addVertexGlobalId(globalId1);
        edge.addVertexGlobalId(globalId2);
        // check edge's existence
        int edgeGlobalId;
        bool nonexist = true;
        for (int j = 0; j < this->vertexGridMeshEdges.size(); ++j)
            if (edge.isDuplicate(this->vertexGridMeshEdges[j])) {
                nonexist = false; // edge exists, use it
                edgeGlobalId = j;
                break;
            }
        if (nonexist) { // edge doesn't exist, add one
            edgeGlobalId = this->vertexGridMeshEdges.size();
            // calculate unit normal vector of the edge here!!!
            this->vertexGridMeshEdges.push_back(edge);
        }
        cell.addEdgeGlobalId(edgeGlobalId);
    }
    this->vertexGridMesh.push_back(cell);
}

template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::
createEdgeCntrGridCell(const CellType cellType,
                     const vector<Grid<ET> >& vertices,
                     int* vertexConn)
{
    Cell cell;
    // set vertices
    for (int i = 0; i < vertices.size(); ++i) {
        // check vertex's existence
        int vertexGlobalId;
        bool nonexist = true;
        for (int j = 0; j < this->edgeCntrGrid.size(); ++j) // This should be optimized
            if (vertices[i].isDuplicate(this->edgeCntrGrid[j])) {
                nonexist = false; // vertex exists, use it
                vertexGlobalId = j;
                break;
            }
        if (nonexist) { // vertex doesn't exsit, add it
            vertexGlobalId = this->edgeCntrGrid.size();
            this->edgeCntrGrid.push_back(vertices[i]);
        }
        cell.addVertexGlobalId(vertexGlobalId);
    }
    // set edges
    int numEdge;
    switch (cellType) {
        case Triangular:
            numEdge = 3;
            break;
        case Quadrilateral:
            numEdge = 4;
            break;
        default:
            break;
    }
    for (int i = 0; i < numEdge; ++i) {
        int localId1, localId2;
        int globalId1, globalId2;
        localId1 = *(vertexConn+i*2);
        localId2 = *(vertexConn+i*2+1);
        globalId1 = cell.getVertexGlobalId(localId1);
        globalId2 = cell.getVertexGlobalId(localId2);
        Edge edge;
        edge.addVertexGlobalId(globalId1);
        edge.addVertexGlobalId(globalId2);
        // check edge's existence
        int edgeGlobalId;
        bool nonexist = true;
        for (int j = 0; j < this->edgeCntrGridMeshEdges.size(); ++j)
            if (edge.isDuplicate(this->edgeCntrGridMeshEdges[j])) {
                nonexist = false; // edge exists, use it
                edgeGlobalId = j;
                break;
            }
        if (nonexist) { // edge doesn't exist, add one
            edgeGlobalId = this->edgeCntrGridMeshEdges.size();
            // calculate unit normal vector of the edge here!!!
            this->edgeCntrGridMeshEdges.push_back(edge);
        }
        cell.addEdgeGlobalId(edgeGlobalId);
    }
    this->edgeCntrGridMesh.push_back(cell);
}

template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::
createCellCntrGridCell(const CellType cellType,
                       const vector<Grid<CT> >& vertices,
                       int* vertexConn)
{
    Cell cell;
    // set vertices
    for (int i = 0; i < vertices.size(); ++i) {
        // check vertex's existence
        int vertexGlobalId;
        bool nonexist = true;
        for (int j = 0; j < this->cellCntrGrid.size(); ++j) // This should be optimized
            if (vertices[i].isDuplicate(this->cellCntrGrid[j])) {
                nonexist = false; // vertex exists, use it
                vertexGlobalId = j;
                break;
            }
        if (nonexist) { // vertex doesn't exsit, add it
            vertexGlobalId = this->cellCntrGrid.size();
            this->cellCntrGrid.push_back(vertices[i]);
        }
        cell.addVertexGlobalId(vertexGlobalId);
    }
    // set edges
    int numEdge;
    switch (cellType) {
        case Triangular:
            numEdge = 3;
            break;
        case Quadrilateral:
            numEdge = 4;
            break;
        default:
            break;
    }
    for (int i = 0; i < numEdge; ++i) {
        int localId1, localId2;
        int globalId1, globalId2;
        localId1 = *(vertexConn+i*2);
        localId2 = *(vertexConn+i*2+1);
        globalId1 = cell.getVertexGlobalId(localId1);
        globalId2 = cell.getVertexGlobalId(localId2);
        Edge edge;
        edge.addVertexGlobalId(globalId1);
        edge.addVertexGlobalId(globalId2);
        // check edge's existence
        int edgeGlobalId;
        bool nonexist = true;
        for (int j = 0; j < this->cellCntrGridMeshEdges.size(); ++j)
            if (edge.isDuplicate(this->cellCntrGridMeshEdges[j])) {
                nonexist = false; // edge exists, use it
                edgeGlobalId = j;
                break;
            }
        if (nonexist) { // edge doesn't exist, add one
            edgeGlobalId = this->cellCntrGridMeshEdges.size();
            // calculate unit normal vector of the edge here!!!
            this->cellCntrGridMeshEdges.push_back(edge);
        }
        cell.addEdgeGlobalId(edgeGlobalId);
    }
    this->cellCntrGridMesh.push_back(cell);
}
#endif

// ************************ mesh interface for external computational components

template <typename VT, typename ET, typename CT>
int MeshManager<VT, ET, CT>::getNumCell(const int dimId) const
{
    if (this->meshStructure != Structured || this->meshForm != Uniform)
        throw runtime_error("MeshManager::getNumCell: "
                            "Invalid for unstructured mesh.");
    if (dimId < 0)
        throw invalid_argument("MeshManager::getNumCell: "
                               "Dimension id must not be less than 0.");
    if (dimId >= this->numDim)
        throw runtime_error("MeshManager::getNumCell: Wrong dimension id.");
    return this->numCell[dimId];
}

template <typename VT, typename ET, typename CT>
int MeshManager<VT, ET, CT>::queryLocation(const vector<double>& x) const
{
    int index = -1;
    if (this->meshForm == Uniform) {
        int o[this->numDim];
        for (int i = 0; i < this->numDim; ++i) {
            o[i] = static_cast<int> (x[i]/this->dx[i]) + 1;
            if (o[i] > this->numCell[i]) {
                o[i] = o[i]-1;
            }
            
        }
        for (int i = this->numDim-1; i > 0; --i) {
            int l = 1;
            for (int j = i-1; j > 0; --j) {
                l *= this->numCell[j];
            }
            index += o[i]*l;
        }
    }
    return index;
}

template <typename VT, typename ET, typename CT>
void MeshManager<VT, ET, CT>::checkBoundary(vector<double>& x) const
{
    if (this->meshStructure == Structured && this->meshForm == Uniform) {
        for (int i = 0; i < this->numDim; ++i) {
            if (this->bndType[i][0] == Periodic) {
                if (x[i] < this->bndX[i][0]) {
                    x[i] += this->span[i];
                }
            }
            if (this->bndType[i][1] == Periodic) {
                if (x[i] > this->bndX[i][1]) {
                    x[i] -= this->span[i];
                }
            }
        }
    }
    return;
}